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Abstract. We consider the leading order quasicontinuum limits of a one- 
dimensional granular medium governed by the Hertz contact law under pre- 
compression. The approximate model which is derived in this limit is justified 
by establishing asymptotic bounds for the error with the help of energy esti- 
mates. The continuum model predicts the development of shock waves, which 
are also studied in the full system with the aid of numerical simulations. We 
also show that existing results concerning the Nonlinear Schrodinger (NLS) 
and Korteweg de-Vries (KdV) approximation of FPU models apply directly to 
a precompressod granular medium in the weakly nonlinear regime. 

1. Introduction. We consider a one-dimensional granular medium which is gov- 
erned by the Hertz law. Denote Qn the relative displacement from equilibrium of 
the n-th particle. Then the renormalized equations of motion describing the q„'s 
have the form [12], 

qn^W'{qn-i-qn)-W'{qn-qn+l), (neZ), (1) 

where 

W'iu) = [So + u]^ , M+ = u e{u) = max(0, u) (2) 

where is the Heaviside function and Sq is the static load (precompression) applied 
to the chain at all times. There is a double nonlinearity stemming from this system 
due to the lack of tensile strength (resulting in an asymmetric potential function 
W) and the nonlinear coupling. For spherical particles we have p = 3/2, which 
corresponds to the classical Hertz contact law. 

It will be convenient to work with the difference u„ = qn-i — q-n, i-C. the strain, 
which satisfies 

Un = W'iun+i)-2W'iu„) + W'iun-i), (neZ). (3) 
With the long wave ansatz 

u„(t) = A{X, T), X = en, T = et 
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where X, T, A{X, T) e M and e ^ 1 is a small perturbation parameter, one can 
derive the leading order continuum model approximation to the discrete dynamics 

d^A^d],{{So+A)P). (4) 
The main result of the present paper is the following approximation theorem. 

Theorem 1.1. Fix Sq > and let peR+ and let A e C ([0, To] , H'^) with 
sup sup |A(X,r)| < (5o/2 and sup || r)||H4 < Ci 

T£[0,To]Xm T£[0,To] 

be a solution of (4). Then for Ci > sufficiently small there exists C,So > such 
that for all e G (0,eo) there are solutions (u„(t))ngz of (3) satisfying 

sup sup \un{t) - A{en, et)\ < Ce^/^. 

te[0,To/e] neZ 

In the weakly nonlinear case, i.e. if the precompression is much greater than 
the amplitude of the solution, then the Nonlinear Schrodinger (NLS) and Korteweg 
de-Vries (KdV) equations can be derived as continuum models for solutions that 
have amplitude of order 0{e) and O(e^) respectively, see Section 8. For solutions 
of order 0(1) (i.e. solutions with amplitude that are on the same order as the 
precompression) model (4) is relevant. In order to establish local existence and 
uniqueness of the continuum model and to handle the nonlincarity in Fourier space, 
we will expand the nonlinearity as a series. Since the amplitude is 0(1) we cannot 
truncate this series (which is done when deriving the NLS and KdV equations). 
Therefore, we require that the precompression is of the same order as the amplitude 
(but not greater). It is unclear if this is only a technical assumption or if it is really 
necessary for having the approximation property. It will be the subject of future 
research to explore whether or not the assumptions of Theorem 1.1 can be relaxed, 
and if so, how the continuum model (4) needs to be modified. 

If one does not ignore higher order terms in the derivation (see next section), 
then one arrives at the following 0(1) continuum model for the strain in the purely 
nonlinear case (if So = 0) 

d'^A = d],iAn + ^dt{An (5) 

which is the equation derived in [1] (where the small parameter e is formally set to 
unity). From a physical standpoint, the reason for doing so is understandable, as 
the higher order term in the model affords the existence of exact localized traveling 
wave solutions (which do not exist in the continuum model (4)). An alternative 
approach towards such a (higher order) quasi-continuum model is the one spear- 
headed in the earlier work of [12], where one does a similar formal Taylor-expansion 
based calculation at the displacement level and then differentiates with respect to x 
to derive the effective long- wavelength equation for the strains r = qx- Remarkably, 
it should be pointed out, as indicated in [1], that these two procedures (reverting 
to strain variables and Taylor expanding to go to long wavelengths) do not com- 
mute, a feature which poses a mathematical challenge in its own right (about their 
respective validity). It should be noted that from a visual inspection the solitary 
wave profiles predicted by (5) (and the corresponding ones of [12]) compare quite 
well with the "numerically exact" traveling wave solution of the granular crystal 
model (3) [1]. Nevertheless, without an approximation theorem, it is not clear to 
what level/extent these kinds of approximations can be used. For example, for the 
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initial value problem (and not just for very special solutions, such as the exact trav- 
eling wave ones), there are several nontrivial concerns about models such as those 
of [1, 12] in connection to their local uniqueness and existence properties and the 
fact the solutions will depend on the small parameter e (hence it is unclear what, 
if any, influence the 0{e) terms have on the dynamics on the 0(1) time scale w.r.t. 
T). Thus, in terms of providing rigorous error estimates for the initial value prob- 
lem, model (4) is more appropriate, at least as a starting point or a leading order 
approximation. Although Eq. (4) fails to describe exact traveling solitary waves 
(and hence the proper, controllable generalization to higher order so as to capture 
this important trait is, from a rigorous perspective, an open problem), the model 
does well in other regards, such as in the description of shock wave formation. This, 
and other properties of (4) are described in Sections 4 and 7. 

Notation. Throughout this paper many possibly different constants arc de- 
noted with the same symbol C if they can be chosen independently of the small 
perturbation parameter < e <C 1. 



2. Fourier transform as fundamental tool. The Fourier transform is the major 
tool in the proof of the approximation result. In this section we recall some basic 
facts and establish notation conventions. 

Since the solutions u of (3) live on Z we need the Fourier transform on Z leading 
to periodic functions in Fourier space. 

Fourier transform on Z: System (3) can be transferred into Fourier space by 

u(fc, t) = J-(7.) (fc, t) = ^ ^ w„ (i)e-'*'-" . (6) 

nGZ 

The inverse of is given by 

u^{t) = {T-^uUt)= [ u{k,t)e^'-dk. (7) 

J —IT 

For every s > the Fourier transform T is continuous from 

£l = {u:Z^R\ \\u{-)\\e2 < oo} 

into 



7Jpgj.(M, C) ~{u : M C I u{-) is s times weakly differentiable, 

\\u{-)\\h^ < oo, u{k) = u{k + 2Tr)} 



where 



= Ed"" ^^^d ii^i(-,t)iii,.^,^ = E( r \dj:'u{k,t)\'dk). 

neZ, m=0 -^^^ 

The inverse Fourier transform J-"^^ is continuous from H^^^ into and from L^^^ 
into i°° where 

= / \u{k,t)\dk and ||u(-, i)||£^ = sup |u,i(t)|. 
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Fourier transform on M: Beside tlie Fourier transform on Z we need the Fourier 
transform on tlic real line in order to handle the continuum model which lives on 
the real line. Wc set 

1 f°° 

u{k,t) ^ I-{u){k,t) ^ — u{x,t)e-'^''dx. (8) 
The inverse is given by 

/oo 
u{k,t)e'^''dk. (9) 
-OO 

For every 5 > the Fourier transform is continuous from 
Ll = {u:R^<C \ <oo} 

into 

i/*(M,C) = {m : M C I u{-) is s times weakly differentiable, < oo } 

where 

{\u{x,t)\^{l+xy)dx and \\u{;t)\\l^ = Y.^ Wu{k,t)?dk). 

The Fourier transform J- and its inverse J-~^ are continuous from to and 
vice versa. Moreover, they are continuous from L]. into C| where 



\u{- 



/OO 
\u{k,t)\{l + ky/^dk and \\u{-,t)\\c^^^^sup\diu{x,t)\ 
-00 j-_Q aSR 



3. Derivation of the continuum model. Wc assume the existence of a > 
such that inf„gzw„(t) > —Sq for all t > 0. This is equivalent to requiring that 
the beads remain in contact for all times. The existence of such solutions will be 
justified below. With this assumption we may ignore the Heaviside function in the 
potential, i.e. we have W'{u) = {So + u)^. Thus, we may expand the nonlinearity 
in a Taylor series, 

OO 

with real- valued coefficients bi = bi{So,p). With this expansion we find 

OO 

= J2 - 2< + (10) 

e=i 

Taking the Fourier-transform of the right hand side of (10) yields 



e ikr, 

u„e 



OO 



= -ujikf E E ^^""6^''" = -^(^)' E ^^^*^(^) 

n6Zf=l e=l 
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where u*^ denotes the {i— l)-tinies convolution of w with itself, ijj{k)'^ = 2(1 — cos(fc)) 
and u{k, t) = u{k + 27r, t). In Fourier space, the system (3) is therefore given by 

oo 

dMk, t) = -uo{kf biu*\k). (11) 

£=1 

We should note here in passing that the use of Fourier space techniques in the 
context of granular systems was pioneered in the work of [3], where it was used to 
develop an understanding of the decay properties of traveling waves in the absence 
of precompression, as well as to offer an efficient numerical tool for computing them. 
Among the recent ramifications of this approach are the proof of the existence of 
such bell-shaped traveling waves without [20] and with [19] precompression. 
The long wave limit ansatz u(n, t) = A(en, et) is given in Fourier space by 

u{k,t)^e-'^A{K,T), = T = et (12) 

with yl : M ^ C a function decaying to zero for — >■ oo. Inserting this ansatz into 
(11), rescaling the integrals, taking formally the limit /^^y^ ^ yields 

oo 

d^A{K, T) = -K^ biA*\K, T) + 0(e2) 

£=1 

where we used that cj(fc)^ = fc^ + 0{k^) = e'^K'^ + 0{{eK)^). Ignoring the higher 
order terms and taking the inverse Fourier transform of this expression yields our 
continuum model (4). Note that keeping the next term in the expansion of cj^ would 
yield Eq. (5) which, as mentioned above, is not covered by the proof presented 
herein. 

Before we turn these formal calculations into rigorous arguments we consider 
properties of the continuum model (4) itself. 

4. Local existence and uniqueness of the continuum equation. Wc will 
need a certain regularity of the solutions of the continuum model (4), therefore, we 
prove the following existence and uniqueness result for the limit equation. 

Lemma 4.1. Fix (5o > and s > 4. Let Aq = ^o(-) G H'' satisfy sup^g^ |Ao(X)| < 
Sq and ||j4o||^s < Ci for a Ci > 0. Then for Ci > sufficiently small there exists 
a Tq > and solutions A ~ A{X, T) of (4) with 

A e C([0,To] ,i7'') and sup sup |^(X, T)| < 5o/2. (13) 

Te[o,To] xeR 

Proof. We know A satisfies 

A = dl ((<5o + AY) . 
As before, we expand the nonlinearity in a series 



with real-valued coefficients hi as defined before. The series is convergent for |A| < 
(5o. Applying 9-^ to (4) then multiplying the resulting equation with 9^ dxA and 
then integrating w.r.t. X yields 
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{d^'drA) {dJ^'d^A) dX = / {d^^'drA) Ox {{So + A)P) dX 

-OO J —OO 

{d],^dTA)Y,bidx{A')dX 

-OO 

which yields 

- / dr ({d^^dTA)^] dX = - {dTA)Y,b,A' dX 

J— OO J —OO /)_-, 



£^1 

OOO OO 



/OO J-^ 1 



and so 

drEo = 



where 



Eo = 

Proceeding similarly for the derivatives yields 



/OO poo 
{d'xdrA) (d^d^A) dX = / {d'xdrA) d'^dl ((<5o + Af) dX 
-OO J —OO 

OO 

{dj,dTA)Y,bid'x^^{A')dX. 



— OO 
OO 



By partial integration we obtain 

-| pOO pOO o^ 

^ / Ot {d^drAf dX = - {d^x^'drA) ^ hd'+\A') dX 

^ J —OO J —OO D_-i 



1^1 



1 pOO ^ 

/ ^£M'-i5T((ai+^A)') dX + Gs 

J~oo 

1 poo I ^ \ 



where G^, Gs only contain terms with at most s + 1 spatial and temporal derivatives. 
Therefore 

drEs+i = Gs 

where 



Es+i = 

and 



i J ^ (dl^dTAf + (Y.^b,A'^' {dx'-'Afj dX 



|G,|<G(P||^,,, + ||aTA||^. 
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if ll^ll^s+i and ||9T^||/fs are smaller than half of the radius of convergence of the 
involved series. Note that all series have the same radius of convergence, namely 
So- 

Since a multiple of £s+i = Eq + ■ ■ ■ + Eg+i is an upper bound of the squared 
if+^-norm we obtain an estimate 

Hence by Gronwall's inequality we can guarantee that stays in between 

half the radius of convergence for all t € [0, Tq] if Tq > and Ci > are chosen 
sufficiently small. Since for a: G M the sup-norm can be estimated by the iJ^-norm, 
the second inequality of (13) follows too. Since the previous a priori estimates 
guarantee that we have a quasilincar system in the sense of [7], the local existence 
and uniqueness of solutions follows. □ 



5. Estimates for the residual. For the proof of the approximation result we need 
a way to measure how ansatz (12) fails to satisfy (3), i.e. we will need estimates for 
the residual. 

It turns out to be advantageous to work in Fourier space, i.e., to work with (11) 
instead of (3). The error e/^R ^u- A with A{k, t) = e~^A{K, T) satisfies 



e^d^Rik, t)=- u{kf ^ bi{A + 7^R)*^ik, t) - d'iA{k, t) 
e=i 

oo oo 

= - uj{k)^ be{A + 7^R)*^{k, t) + oj{kf beA*^{k, t) 
1=1 1=1 

oo oo 

;(fc)2 ^ biA*\k, t) + eY, beA*^{k, t) 



and so 



e^d^R{k,t) ^ -iu{kfYbi{A+7PR)*'^{k,t)+w{kfYbeA*\k,t)+R^){k,t) 

i=i i=i 

(14) 



where 



R^^)(yfc, t) = (fc2 - uj{kf) J2 beA*\k, t) 

1=1 

stands for the residual terms, i.e., for the terms which we neglected within the 
leading order quasicontinuum approximation. In order to solve (14) with boundary 
conditions i?(fc, t) = R{k + 2-k, t) we have to modify A{k, t) which decays to zero 
for \k\ — >■ oo. We multiply A{k,t) with a cut-off function X[-7r/2,7r/2]- The segment 
from [— TT, tt] is then extended periodically with a period 27r to the entire real axis. 
Call the outcome .A#(fc, t). In exactly the same way we modify — fc^ in the residual, 
which is then denoted with —k\. 
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We will need an estimate of the difference A{k, t) -~A^{k, t) and the error e^R = 
u — A^ which satisfies 

CJO OQ 

e'^d'iR{k,t) = -uj{kf^bi{AiPr^R)*^{k,t)^^u{kf^hiA*l{k,t)+R^^^#){k,t) 
1=1 1=1 * 

where 

oo 

RiK^)(fc,i) - {kl-u{kf)Y,b,A*l{k,t). 
* 1=1 
Thus, in order to bound the error i?, we will need estimates for the residual 

Res(^^). Since 1? is closed under convolution on [— tt, tt] it turns out to be sufficient 

# 

to make the estimates in I? . We have 

Lemma 5.1. Fja;(5o>0. Lef A G C ([0, Tq] , iJ'*) wit/i supj,g[(,.To] suPxgr l^l-'^) < 
(5o/2 he a solution (4). Then there exist so^C > such that for all e S (0,eo) we 
have 



sup 

te[0,To/s] 



R|s(^#) 



and 



sup 

te[0,To/e] 



UJ ^Res{A#] 



< 



Proof. For completeness wc recall the proof of [2, Lemma 3.3]. We have that Q = 
^eL2biA^{k,t) satisfies 



sup 



1/2 



\g(k,t)\''{i+'-rdk] 



for e 0, (note the index of Q starts at ^ 2). The loss of e~^^^ comes from the 
scaling properties of the L^-norm. Using that k^ — Ld{k)'^ = 0{k'^) then yields 



< C sup 



1/2 

\{k^ -uj{ky)g{k,t)\'-'dk] <c[i \k'^g{k,t)\^dk 
k* 



— TT 

2 \ 1/2 



1/2 



\g{k,t)\-'{i + —)^dk] <Ce^/2 



Since bi is independent of e and since fc^ — uj{k)'^ ~ 0{k^), it follows that 



{e -Lj{k)^)J2btA*jik,t) 



2 \ 1/2 

dk I < Ce^/^ 



Applying the same argument as above with (fc^ — u{kY) / uj{k) = 0{k^) yields 

oo 

-\k){k^ - u{kf)Y,b£A*l{k,t) 



2 \ 1/2 

dk\ <Ce^/2_ 



□ 



We close this section with an estimate for the difference A{k,t) — A^{k,t). 
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Lemma 5.2. Let A e C{\0,To\,H^) then 



sup sup \A{n,t) ~ A#{n,t)\ < 

te[0,To/£] nGZ 



Proof. Wc have 



A{n,t) 



e-^A{-,et)e"'"dk and A#in,t) 

1 ^ 



-7r/2 ^ 



For the difference we obtain 



it/2 



< C sup 

fee[7r/2,cx)] 



(1 + 



\2 



y |£-M(^,rf)|2(l + ^)4dfcj <C£^/2 



uniformly in n. The loss of e again comes from the scaling properties of the 
L^-norm. 

□ 

6. The error estimates. It remains to bound the solutions of (15). In accordance 
with Lemma 5.1 we choose (3 = 3/2. We proceed as in Section 4 using energy 
estimates. 

From Section 5 we know 



*£ 



^ b,{[A# + e^R) -A# -)\ e-P + e-^Res(^#). 



Define the energy 

/TT /"TT 
\uj-^dtR\'^dk/2+ / \uj-^hiR\^dk/2. 
-TT J — TT 

Making use of tlie fact 

Re / uj-^dtRd^Rdk ^ dt \uj-^dtR\^dk/2. 



we can compute 

dtEo^Re\-e-P f idM) \ Y,h{(A#+ePR)'~A#* 



\l=2 



dk 



LJ-'dtR] ( uj-'ResiA#) ] dk. 



Note that the autonomous linear terms have canceled, explaining why the sum 
begins at £ = 2. Recall that the application of o;""'^ to the residual terms is well 
defined, see Lemma 5.1. We show below (see e.g. (15)) that the application of 
Lj^^ to dtR and R is also well defined. Using the Plancherel's identity allows us to 
rewrite 

£-^Re J {dtR)\Y2^bei(A^ + ePRy -A^*^)jdk 
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as 



n& 1=2 j=l 



ee^^eC) -^5. 



^^EIE^^eO-^t^^"'"')^""^^'""^ 



oo t 



l\ 1 1 + 1 N ^ / -/--IN ,3(j-l) 



EE^^E -tt(^"^"')^*(<"^)^' 



where ^„ (resp. i?„) is the inverse discrete Fourier transform of (resp. R) 
evaluated at n. This motivates the definition of a modified energy 

which by construction satisfies, 



dtEi =Ho + e-" J (^^-^dtR) (^lj-^R^^#)^ 

Ho = EE^^E ( ■) -^^iRn^^')OMi-ne^^^ 

neZ£=2 j = l ^ 



Let \\A{-,T)\\p < Ca and \\dtA{-,T)\\p < Cse and \\R{-,t)\\p < Ce with Cb defined 
below in (15). Then 



< {e\\R\\(2Ci+e"+^\\R\\%C5{CE)) 

where 



oo _j 

^4:=E^^2^^~^^^2"'C'3, 
C,{Ce) :=E&.E (!■) -^^C=E\l-^)C'2-'-'C,ef^''-''- 

1=2 j=2 ■' 

The first series is convergent for Ci (and hence C2) sufficiently small. By chang- 
ing indices, factorizing and ignoring a finite number of terms in the second series 
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results in the expression C^{Ce) ~ Sq{1 + (C2 + e^C£;)/((5o))^ and thus the series 
is convergent if Ci and e arc sufficiently smaU. Due to the Plancherel identity the 
energy Ei is an upper bound for the squared L^-norm for Ci > sufficiently small. 
Then since lo is bounded and since C on bounded domains (as a consequence 
is closed under convolution), we find 



Ho < {C4 + l)E,s 
where we pick e > small enough such that e^C^{CE)CE < 1- Thus 

dtEi < (C4 + 1 + C(i)E^e + Cee 
with constants Cj independent of s. For the inequality above we used the fact that 

J (^uj-^dtR) (^uj-^R^{A#)^ dk < Cee (^^ + J \^^^^dtR\^dk^ , 

which is a consequence of the Cauchy-Schwarz inequality, Lemma 5.1 and the fact 
that ^/x < 1 + X. Hence by Gronwall's inequality we have 

-P E,it) < ,ic.^^^c.)T. _ d (15) 

tg[0,To/e] 04 + 1 + Og O4 + i + Og 

for aU t e [0, To/e] where we set u(-, 0) = A{e-,0) such that i?i(0) = 0. As a direct 
consequence we have 

sup sup|M„(t) -yl#(n,i)| < Ce^/^. 

te[0,To/e] nSZ 

Combining this estimate with the one from Lemma 5.2 gives the assertion of The- 
orem 1.1, completing our proof. 

□ 

7. Shock Formation in Granular Media. We now briefly discuss some analyt- 
ical and numerical observations that stem from the leading order quasicontinuum 
approximation developed in Theorem 1.1. The continuum model 

dlA = dl{{5o + AY) 
can be written as a system of conservation laws 



(16) 



OtA — dxv = 

drv - dx[{5Q + AY] =0 

These conservation laws have the form of a so-called p-system [10, 18]. Defining 
we can express (16) as, 

drU + dxF{U) = 

or equivalently as 

drU + DF{U)dxU ^0 (17) 

where 
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Figure 1. Development of a shock-like structure with a smooth 
and localized initial condition of the granular crystal model (3) 
with p = 3/2 and S ~ 0.1. A space-time contour plot is shown on 
the left and spatial profiles at time t = (solid black line), t = 370 
(red dashed line) and t ~ 750 (blue dashed-dot line) are shown 
in the right panel. In both panels it can be seen that the pulse 
separates into two counter-propagating waves. A zoom of the wave 
before and after the point of shock development is shown in Fig. 2. 



which has the eigenvalues X±{U) — ±y^p{So + A)p^^. Thus, solutions of (4) will 
consist of two counter-propagating waves traveling with velocity ± y/p{So'+~A)P^ . 
Since the wave speed will depend on the amplitude of the solution, bell-shaped initial 
data will deform and steepen and a shock wave will form in finite time. Discrete 
shock-like structures (which we will simply call shock waves) have been studied 
in granular media, e.g. in homogeneous and periodic chains with p = 3/2 and 
(5o = [5, 11]. In those works, however, the shock wave is generated by applying 
a velocity to a single bead [5] or by imparting velocity to the end of the chain 
continuously [11]. In this paper, the mechanism for the development of the shock 
wave is fundamentally different. It manifests from an arbitrary non-monotonically- 
increasing initial strain profile under precompression and given a sufficiently long 
time to develop. 

A natural question is if the shock wave formation predicted by the continuum 
model (4) is also present in full system (3). Theorem 1.1 no longer applies in this 
case, as the the shock wave violates the required smoothness condition. Nonetheless, 
we carry out a numerical simulation to address the relevant question with a smooth 
initial condition. Figure 1 shows the development of a shock wave in the discrete 
model (3) with the initial condition Un(0) = a scch(5en), Vn{0) ~ where a, & G M are 
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t= 245 t= 545 t= 1 000 




Figure 2. A zoom of the wave before (left) and at the approxi- 
mate moment of (middle) and after (right) the point of shock de- 
velopment. The dark solid line is from the direct simulation of the 
discrete model, and the dashed red line is the prediction based on 
the continuum model. Notice that past the point of expected wave 
breaking, the lattice model develops microscopic oscillations. The 
profile at t = 145 was used to construct the red-dashed line using 
the velocity relationship \Jp{5o + A)p~^ , see text. Before the de- 
velopment of the shock, this prediction is very accurate, see e.g. 
the left panel. 

shape parameters. The wave propagation closely follows the theoretical expectation 
on the basis of the leading order quasicontinuum approximation. Moreover, the 
predicted velocity ±-\/p((5o + A)p^^ proved to be very accurate (for example, there 
was less than a %0.01 relative error in the case shown in Fig. 1). 

Wc would like to make a direct comparison of the solutions of the continuum 
model and the discrete model. However, one has to be careful when using nu- 
merical approximations of the continuum model. For example, if one uses a finite 
difference approximation for a spatial discretization of (4), then we arrive at a 
model identical to the granular model (3). One could use other numerical schemes, 
such as those based on adding artificial dispersion [10], but we will proceed in an 
alternative way to predict the development of a shock wave. Using the velocity 
relationship ±\/p(5o + A)p~'^ , we can construct the profile of the continuum model 
for an arbitrary time, starting with the left or right wave (once they are separated) 
as an initial profile, see Fig. 2. A numerical computation is used to separate the 
profiles (i.e. we simulate (4) with a finite difference method until separation but 
before the development of any shock wave, thus avoiding any issues with smooth- 
ness). The continuum model predicts a shock wave for any time past the point of 
non-single- valuedness, as shown in right panel Fig. 2. 

In FPU lattices, it is well known that dispersive shocks can develop, in which 
microscopic oscillations spread out in space and time [6]. From Fig. 2 one clearly 
sees near the point of wave breaking the development of such oscillations, which are 
absent in the continuum model [6]. Thus, it would be relevant to extend works like 
[6] in order to better understand shock waves in the granular crystal model (3). 

8. Justification of the KdV and NLS approximation. In this section we 
would like to contrast the previous result with approximation results for the KdV 
and NLS approximation. The major difference lies in the ratio between the ampli- 
tude and the precompression. For the KdV and NLS approximation this ratio is 
©(e^) resp. 0{e) where < e ^ 1 is the small perturbation parameter, whereas for 
the quasicontinuum approximation the ratio is 0(1), i.e., of a comparable order. 



14 



CHRISTOPER CHONG, P.G. KEVREKIDIS AND GUIDO SCHNEIDER 



Since Eq. (10) is exactly of the form of the FPU systems considered in [2, 15] the 
approximation theorems of these papers also apply here. The subsequent solutions 
Un will be 0{e'^), resp. 0{e) and hence will always live in the ball of the convergence 
with radius 6q if the perturbation parameter < e ^ 1 is chosen to be sufficiently 
small. 

Even spatially periodic arrangements can be considered here, namely 

dfun a„+i«_^i) - 2a„«) + a„_i«_i) (18) 

with a„ = ttn+N for a fixed N. 

We formulate the relevant approximation theorems in the homogeneous case 
= 1. For the KdV approximation we have, 

Theorem 8.1. Let A E C([0,ro] ,H^) be a solution of the KdV equation drA ~ 
vid\A + V2dx{A^) with suitable chosen coefficients i'i,V2 € K- Then there exist 
£0 > 0, C > such that for all e £ (0,eo) we have solutions {un)n£i of (18) with 

sup sup \un{t) - V'„(i)l < C'e^/^ 

te[0,To/e3] riGN 

where 

Mt) =e'^A{e{n-uj[{0)t),eh) 

with u}i{k)'^ = uj{kybi. 

Proof. The proof follows trivially from [2, Theorem 3.1] or [16]. We note there 
is no gap-opening in the small amplitude limit (i.e. the Heaviside function in the 
potential will play no role). □ 

The theorem can be generalized easily to an approximation theorem for two 
decoupled KdV equations describing counter-propagating waves, cf. [16]. 
For the NLS approximation we have, 

Theorem 8.2. Let A G C{[0,To], H^^) be a solution of the NLS equation drA = 
ividxA-\-iv2A\A^ with suitable chosen coefficients vi,V2 G M. Then there exist 
Eo > 0, C > such that for all e S (0,eo) we have solutions {un)nei. of (18) with 

sup sup \Unit) - ^n{i)\ < CC'I'^ 
te[0,To/e2] n&\ 

where 

V'„(t) = eA {e{n - uj[iko)t) , eh) e'C^"""""*) + c.c. 
with a;i(fc)^ = U!{k)^bi 

Proof. The proof follows trivially from [2, Theorem 4.1] or [15, Theorem 1.1]. □ 

9. Conclusions and Future Challenges. The main result of this paper derives in 
a rigorous way (and with controllable corrections) the leading order quasicontinuum 
approximation for long wavelength solutions in the granular crystal model (3), in 
accordance with the formal derivation of Nestercnko [12] . As a technical assumption, 
we required the presence of a precompression factor (while the original Nesterenko 
model has been developed also in the case of the so-called "sonic vacuum" i.e., 
without precompression). One obvious avenue of future research is to investigate 
a proof without this assumption, which, however, would demand a fundamentally 
different technique than the one presented herein. 

On the other hand, perhaps an even more important aspect of investigation 
concerns the well-posedness theory of the full Nesterenko model [12] or of the variant 
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developed by Ahnert and Pikovsky in [1]. One important consequence of keeping 
only the leading order terms in the continuum model (as done herein) is the inability 
to capture the exact solitary wave solutions (which are known to exist in granular 
crystals [4, 3]) and are at the core of experimental observations in such systems [12, 
17, 8]. The methods of this paper cannot be directly applied to that case, although 
we should note that we suspect that these higher order long wavelength models suffer 
(especially so in the case of precomprcssion) from the type of pathologies that were 
identified by Rosenau and led him to devise appropriate regularizations [13, 14]; see 
also the more recent discussion of [9] . It would be especially relevant to consider such 
regularizations of the higher order long wavelength models both from a rigorous, as 
well as from a numerical perspective. 

Another important consequence of keeping only first order terms in the contin- 
uum model is the prediction of the development of shock waves for a suitable (yet 
broad) class of initial data. Although the main theorem of this paper docs not apply 
to the case of shock waves, due to smoothness considerations, numerical simulations 
indicate the steepening of relevant initial data towards a shock structure and sug- 
gest that this is indeed an issue worthy of further exploration, with an aim towards 
transferring these results to the discrete model. Indeed, it is known in FPU lattices 
that this procedure fails [6] , due to the existence of high frequency oscillations (re- 
sulting from so-called dispersion shock waves). Thus, a different continuum model 
(than the one derived herein) will most likely be needed to fully characterize the 
emerging dispersive shock wave case. These topics are currently under consideration 
and will be reported in future publications. 
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